{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# COVID-19 Simulation\n",
    "\n",
    "This COVID-19 simulation was developed by [transentis consulting](https://www.transentis.com) from Berlin, based on the [SIR model](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model).\n",
    "\n",
    "Please read the companion blog post [Covid 19 and the SIR Model](https://www.transentis.com/covid-19-sir-model/) for details on the model and a discussion of COVID-19 scenarios.\n",
    "\n",
    "The simulation was built using the [BPTK-Py](https://bptk.transentis-labs.com) framework, the dashboard using [Voila](https://voila.readthedocs.io/en/stable/). You can find the complete source code on [GitHub](https://github.com/transentis/sim-covid-19). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "document.title=\"transentis | Covid-19 simulation\""
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name=\"transentis | Covid-19 simulation\"\n",
    "from IPython.display import display,Javascript\n",
    "Javascript('document.title=\"{}\"'.format(name))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "37038eeeb1f94f50baf280b06f0794fa",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Tab(children=(Output(), Output(), Output(), Output(), Output()), _titles={'0': 'Population', '1': 'Intensive C…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8afb40247e39451298229fd23b3b4bab",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "GridBox(children=(IntSlider(value=0, continuous_update=False, description='Visualization Period Begin', layout…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from BPTK_Py import Model\n",
    "from BPTK_Py import sd_functions as sd\n",
    "\n",
    "model = Model(starttime=1.0,stoptime=1500.0,dt=1.0,name='sir')\n",
    "\n",
    "susceptible = model.stock(\"susceptible\")  # those that have not been infected yet\n",
    "infectious = model.stock(\"infectious\") # those that are currently carrying an infection\n",
    "recovered = model.stock(\"recovered\") # those that have recovered from an infection\n",
    "deceased = model.stock(\"deceased\") # those that have died from the infection\n",
    "\n",
    "infection_rate = model.flow(\"infection_rate\") # the rate at which people are becoming infected\n",
    "recovery_rate = model.flow(\"recovery_rate\") # the rate at which people are recovering from an infection\n",
    "death_rate = model.flow(\"death_rate\") # the rate at which people are dying from an infection\n",
    "\n",
    "infectivity = model.constant(\"infectivity\") # the infectivity of the corona virus\n",
    "lethality = model.constant(\"lethality\") # the lethality of the corona virus\n",
    "normal_contact_rate=model.constant(\"normal_contact_rate\") # the average contact rate between people in \"normal\" times\n",
    "duration = model.constant(\"duration\") # the average time it takes to recover from the virus\n",
    "\n",
    "intensive_available = model.constant(\"intensive_available\") # the number of intensive care units available\n",
    "intensive_percentage = model.constant(\"intensive_percentage\") # the fraction of people needing intensive care\n",
    "\n",
    "dashboard_on=model.constant(\"dashboard_on\") # should be equal to 1.0 for dashboard scenario and 0.0 otherwise\n",
    "distancing_on=model.constant(\"distancing_on\") # defines whether social distancing is on within an interactive scenario\n",
    "distancing_contact_rate=model.constant(\"distancing_contact_rate\") # interactive scenario: the contact rate when practicing social distancing\n",
    "distancing_begin=model.constant(\"distancing_begin\") # interactive scenario: when to begin social distancing\n",
    "distancing_duration=model.constant(\"distancing_duration\")# interactive scenario: when social distancing ends and we return to normal behavior\n",
    "\n",
    "total_population = model.converter(\"total_population\") # the total population, i.e the sum of susceptible, infected and recovered\n",
    "contact_rate = model.converter(\"contact_rate\") # the rate at which people are being contacted, in all scenarios\n",
    "\n",
    "intensive_needed = model.converter(\"intensive_needed\") # the number of intensive care units needed at any time\n",
    "\n",
    "contact_number = model.converter(\"contact_number\") #  measures which fraction of a susceptible population is infected by a contagious person\n",
    "reproduction_rate = model.converter(\"reproduction_rate\") # measures the rate at which an epidemic reproduces\n",
    "\n",
    "variable_contact_rate=model.converter(\"variable_contact_rate\") # the variable contact rate used in the detailed scenarios\n",
    "dashboard_with_distancing_contact_rate=model.converter(\"dashboard_with_distancing_contact_rate\") # the contact rate used in interactive scenarios with distancing on\n",
    "dashboard_contact_rate=model.converter(\"dashboard_contact_rate\") # the contact rate used in interactive scenarios with distancing off\n",
    "\n",
    "susceptible.initial_value = 80000000.0\n",
    "infectious.initial_value = 120.0\n",
    "recovered.initial_value = 0.0\n",
    "deceased.initial_value = 0.0\n",
    "\n",
    "infectivity.equation = 0.02\n",
    "duration.equation = 20.0\n",
    "lethality.equation = 0.001\n",
    "\n",
    "intensive_percentage.equation = 0.002\n",
    "intensive_available.equation = 30000.0\n",
    "intensive_needed.equation =infectious*intensive_percentage\n",
    "\n",
    "susceptible.equation = -infection_rate\n",
    "infectious.equation = infection_rate - recovery_rate - death_rate\n",
    "recovered.equation = recovery_rate\n",
    "deceased.equation = death_rate\n",
    "\n",
    "total_population.equation = susceptible+infectious+recovered\n",
    "\n",
    "infection_rate.equation = (contact_rate*infectivity*infectious)*(susceptible/total_population)\n",
    "\n",
    "recovery_rate.equation = infectious/duration\n",
    "\n",
    "death_rate.equation = infectious*lethality\n",
    "\n",
    "contact_number.equation=contact_rate*infectivity*duration\n",
    "reproduction_rate.equation=contact_number*(susceptible/total_population)\n",
    "\n",
    "contact_rate.equation=20\n",
    "\n",
    "contact_rate.equation=dashboard_on*dashboard_contact_rate+(-dashboard_on+1.0)*variable_contact_rate\n",
    "dashboard_on.equation=1.0\n",
    "distancing_on.equation=0.0\n",
    "dashboard_contact_rate.equation=distancing_on*dashboard_with_distancing_contact_rate+(-distancing_on+1.0)*normal_contact_rate\n",
    "dashboard_with_distancing_contact_rate.equation=sd.If(sd.Or(sd.time()<distancing_begin,sd.time()>distancing_begin+distancing_duration),normal_contact_rate,distancing_contact_rate)\n",
    "normal_contact_rate.equation=20.0\n",
    "distancing_contact_rate.equation=2.0\n",
    "distancing_begin.equation=20.0\n",
    "distancing_duration.equation=200.0\n",
    "\n",
    "variable_contact_rate.equation=sd.lookup(sd.time(),\"variable_contact_rate\")\n",
    "\n",
    "variable_contact_rate_points = [[0,20.0],[1500,20.0]]\n",
    "\n",
    "model.points[\"variable_contact_rate\"]=variable_contact_rate_points\n",
    "\n",
    "import BPTK_Py\n",
    "import pandas as pd\n",
    "bptk = BPTK_Py.bptk()\n",
    "bptk.register_model(model)\n",
    "\n",
    "bptk.register_scenarios(\n",
    "    scenarios ={\n",
    "        \"base\": {},\n",
    "        \"weakSocialDistancing\": {},\n",
    "        \"strongSocialDistancing\": {},\n",
    "        \"shortTermMeasures\": {},\n",
    "        \"dashboard\":{}\n",
    "    },\n",
    "    scenario_manager=\"smSir\")\n",
    "\n",
    "bptk.reset_scenario_cache(scenario_manager=\"smSir\", scenario=\"dashboard\")\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import HTML\n",
    "out_population = widgets.Output()\n",
    "out_intensive_care = widgets.Output()\n",
    "out_indicators = widgets.Output()\n",
    "out_contact_rate = widgets.Output()\n",
    "out_assumptions = widgets.Output()\n",
    "\n",
    "wdg_style_distancing_on = {'description_width': 'initial'}\n",
    "wdg_distancing_on = widgets.Checkbox(\n",
    "        value=False,\n",
    "        continuous_update=False,\n",
    "        disabled=False,\n",
    "        description='Introduce Social Distancing',\n",
    "        style=wdg_style_distancing_on\n",
    "    )\n",
    "wdg_distancing_on.layout.width='450px'\n",
    "\n",
    "wdg_style_normal_contact_rate = {'description_width': 'initial'}\n",
    "wdg_normal_contact_rate=widgets.FloatSlider(\n",
    "        value=20,\n",
    "        min=0,\n",
    "        max=40,\n",
    "        step=1,\n",
    "        continuous_update=False,\n",
    "        description='Normal Contact Rate',\n",
    "        style=wdg_style_normal_contact_rate\n",
    "    )\n",
    "wdg_normal_contact_rate.layout.width='450px'\n",
    "\n",
    "wdg_style_vis_begin= {'description_width': 'initial'}\n",
    "wdg_vis_begin=widgets.IntSlider(\n",
    "        value=0,\n",
    "        min=0,\n",
    "        max=1499,\n",
    "        continuous_update=False,\n",
    "        description='Visualization Period Begin',\n",
    "        style=wdg_style_vis_begin\n",
    "    )\n",
    "wdg_vis_begin.layout.width='450px'\n",
    "\n",
    "wdg_style_vis_end = {'description_width': 'initial'}\n",
    "wdg_vis_end=widgets.IntSlider(\n",
    "        value=1500,\n",
    "        min=1,\n",
    "        max=1500,\n",
    "        continuous_update=False,\n",
    "        description='Visualization Period End',\n",
    "        style=wdg_style_vis_end\n",
    "    )\n",
    "wdg_vis_end.layout.width='450px'\n",
    "\n",
    "wdg_style_distancing_contact_rate = {'description_width': 'initial'}\n",
    "wdg_distancing_contact_rate=widgets.FloatSlider(\n",
    "        value=2,\n",
    "        min=0,\n",
    "        max=20,\n",
    "        step=1,\n",
    "        continuous_update=False,\n",
    "        description='Distancing Contact Rate',\n",
    "        style=wdg_style_distancing_contact_rate\n",
    "    )\n",
    "wdg_distancing_contact_rate.layout.width='450px'\n",
    "\n",
    "\n",
    "wdg_style_distancing_begin= {'description_width': 'initial'}\n",
    "wdg_distancing_begin=widgets.FloatSlider(\n",
    "        value=20,\n",
    "        min=0,\n",
    "        max=200,\n",
    "        step=1,\n",
    "        continuous_update=False,\n",
    "        description='Distancing Begin',\n",
    "        style=wdg_style_distancing_begin\n",
    "    )\n",
    "wdg_distancing_begin.layout.width='450px'\n",
    "\n",
    "wdg_style_distancing_duration= {'description_width': 'initial'}\n",
    "wdg_distancing_duration=widgets.FloatSlider(\n",
    "        value=100,\n",
    "        min=0,\n",
    "        max=1300,\n",
    "        step=1,\n",
    "        continuous_update=False,\n",
    "        description='Distancing Duration',\n",
    "        style=wdg_style_distancing_duration\n",
    "    )\n",
    "wdg_distancing_duration.layout.width='450px'\n",
    "\n",
    "def event_handler_distancing_on(change):\n",
    "    update_graphs(change.new, wdg_normal_contact_rate.value, wdg_vis_begin.value, wdg_vis_end.value, wdg_distancing_contact_rate.value, wdg_distancing_begin.value, wdg_distancing_duration.value)\n",
    "    \n",
    "def event_handler_normal_contact_rate(change):\n",
    "    update_graphs(wdg_distancing_on.value, change.new, wdg_vis_begin.value, wdg_vis_end.value, wdg_distancing_contact_rate.value, wdg_distancing_begin.value,wdg_distancing_duration.value)\n",
    "    \n",
    "def event_handler_vis_begin(change):\n",
    "    update_graphs(wdg_distancing_on.value, wdg_normal_contact_rate.value,change.new, wdg_vis_end.value, wdg_distancing_contact_rate.value, wdg_distancing_begin.value,wdg_distancing_duration.value)\n",
    "    \n",
    "def event_handler_vis_end(change):\n",
    "    update_graphs(wdg_distancing_on.value, wdg_normal_contact_rate.value, wdg_vis_begin.value, change.new, wdg_distancing_contact_rate.value, wdg_distancing_begin.value,wdg_distancing_duration.value)\n",
    "     \n",
    "def event_handler_distancing_contact_rate(change):\n",
    "    update_graphs(wdg_distancing_on.value, wdg_normal_contact_rate.value, wdg_vis_begin.value, wdg_vis_end.value, change.new, wdg_distancing_begin.value, wdg_distancing_duration.value)\n",
    "    \n",
    "def event_handler_distancing_begin(change):\n",
    "    update_graphs(wdg_distancing_on.value, wdg_normal_contact_rate.value, wdg_vis_begin.value, wdg_vis_end.value, wdg_distancing_contact_rate.value,change.new,wdg_distancing_duration.value)\n",
    "    \n",
    "def event_handler_distancing_duration(change):\n",
    "    update_graphs(wdg_distancing_on.value, wdg_normal_contact_rate.value, wdg_vis_begin.value, wdg_vis_end.value, wdg_distancing_contact_rate.value,wdg_distancing_begin.value,change.new)\n",
    "   \n",
    "    \n",
    "def update_graphs(distancing_on, normal_contact_rate, vis_begin, vis_end, distancing_contact_rate,distancing_begin, distancing_duration):\n",
    "    scenario= bptk.get_scenario(\"smSir\",\"dashboard\")\n",
    "\n",
    "    scenario.constants[\"distancing_on\"]=1.0 if distancing_on else 0.0\n",
    "    scenario.constants[\"normal_contact_rate\"]=normal_contact_rate\n",
    "    scenario.constants[\"distancing_contact_rate\"]=distancing_contact_rate\n",
    "    scenario.constants[\"distancing_begin\"]=distancing_begin\n",
    "    scenario.constants[\"distancing_duration\"]=distancing_duration\n",
    "    scenario.constants[\"dashboard_on\"]=1.0\n",
    "    \n",
    "    wdg_distancing_contact_rate.layout.visibility = 'visible' if distancing_on else 'hidden'\n",
    "    wdg_distancing_begin.layout.visibility = 'visible' if distancing_on else 'hidden'\n",
    "    wdg_distancing_duration.layout.visibility = 'visible' if distancing_on else 'hidden'\n",
    "    \n",
    "    bptk.reset_scenario_cache(scenario_manager=\"smSir\", scenario=\"dashboard\")\n",
    "\n",
    "    out_population.clear_output(wait=True)\n",
    "    out_intensive_care.clear_output(wait=True)\n",
    "    out_contact_rate.clear_output(wait=True)\n",
    "    out_assumptions.clear_output(wait=True)\n",
    "    out_indicators.clear_output(wait=True)\n",
    "    \n",
    "    with out_population:\n",
    "        # turn of pyplot's interactive mode to ensure the plot is not created directly\n",
    "        plt.ioff()\n",
    "        # clear the widgets output ... otherwise we will end up with a long list of plots, one for each change of settings\n",
    "        \n",
    "        # create the plot, but don't show it yet\n",
    "        plot=bptk.plot_scenarios(\n",
    "            scenario_managers=[\"smSir\"],\n",
    "            scenarios=[\"dashboard\"],\n",
    "            title=\"Recovered Population vs. Deaths\",\n",
    "            x_label=\"Days\",\n",
    "            y_label=\"Persons\",\n",
    "            equations=[\n",
    "                \"infectious\",\n",
    "                \"recovered\",\n",
    "                \"deceased\"\n",
    "            ],\n",
    "             series_names={\n",
    "                 \"smSir_dashboard_infectious\":\"Infectious Population\",\n",
    "                 \"smSir_dashboard_recovered\":\"Recovered Population\",\n",
    "                 \"smSir_dashboard_deceased\":\"Deceased\"     \n",
    "            },\n",
    "            visualize_from_period=vis_begin,\n",
    "            visualize_to_period=vis_end\n",
    "        )\n",
    "        # show the plot\n",
    "        plt.show(plot)\n",
    "        # turn interactive mode on again\n",
    "        plt.ion()  \n",
    "              \n",
    "    with out_intensive_care:\n",
    "        plt.ioff()\n",
    "        plot=bptk.plot_scenarios(\n",
    "            scenario_managers=[\"smSir\"],\n",
    "            scenarios=[\"dashboard\"],\n",
    "            title=\"Available Intensive Care vs. Needed Intensive Care\",\n",
    "            x_label=\"Days\",\n",
    "            y_label=\"Intensive Care Units\",\n",
    "            equations=[\n",
    "                \"intensive_needed\",\n",
    "                \"intensive_available\"\n",
    "            ],\n",
    "            series_names={\n",
    "                \"smSir_dashboard_intensive_needed\": \"Intensive Care Needed\",\n",
    "                \"smSir_dashboard_intensive_available\": \"Intensive Care Available\"  \n",
    "            },\n",
    "            visualize_from_period=vis_begin,\n",
    "            visualize_to_period=vis_end        \n",
    "        )\n",
    "        plt.show(plot)\n",
    "        plt.ion()\n",
    "        \n",
    "    with out_indicators:\n",
    "        plt.ioff()\n",
    "        plot=bptk.plot_scenarios(\n",
    "            scenario_managers=[\"smSir\"],\n",
    "            scenarios=[\"dashboard\"],\n",
    "            title=\"Indicators\",\n",
    "            x_label=\"Days\",\n",
    "            y_label=\"Measure\",\n",
    "            equations=[\n",
    "                \"contact_number\",\n",
    "                \"reproduction_rate\"\n",
    "            ],\n",
    "            series_names={\n",
    "                \"smSir_dashboard_contact_number\": \"Contact Number\",\n",
    "                \"smSir_dashboard_reproduction_rate\": \"Reproduction Rate\" \n",
    "            },\n",
    "            visualize_from_period=vis_begin,\n",
    "            visualize_to_period=vis_end\n",
    "        )\n",
    "        plt.show(plot)\n",
    "        plt.ion()\n",
    "    with out_contact_rate:\n",
    "        plt.ioff()\n",
    "        plot=bptk.plot_scenarios(\n",
    "            scenario_managers=[\"smSir\"],\n",
    "            scenarios=[\"dashboard\"],\n",
    "            title=\"Contact Rate\",\n",
    "            x_label=\"Days\",\n",
    "            y_label=\"Average People Contacted\",\n",
    "            equations=[\n",
    "                \"contact_rate\"\n",
    "            ],\n",
    "            series_names={\n",
    "                \"smSir_dashboard_contact_rate\": \"Contact Rate\" \n",
    "            },\n",
    "            visualize_from_period=vis_begin,\n",
    "            visualize_to_period=vis_end\n",
    "        )\n",
    "        plt.show(plot)\n",
    "        plt.ion()\n",
    "    with out_assumptions:\n",
    "        display(HTML(\"<p>The implementation here is roughly calibrated to the current situation in Germany (as of 27.3.2020). It illustrates the effects of social distancing in achieving the objective of keeping the strain on the health care system as small as possible.</p><ul><li>Contact Rate: 20 persons. Defines how many people a person encounters per day in average.</li><li>Infectivity: 2%. Defines the probability that a person becomes infected after contact with an infectious person.</li><li>Duration. Defines how long an infective person remains contagious</li><li>Population. The susceptible population starts at 80 Mio., the infectious population starts at 120 persons.</li><li>Intensive Care Needed: 0.2%. Measures the number of infected people who need intensive care.</li><li>Intensive Care Available: 30,000 units. The number of intensive care units available.</li></ul><p>With the above settings, this means we have a <i>contact number</i> of 8 in the base settings. The contact number is the product of contact rate, infectivity and duration.</p>\"))\n",
    "        \n",
    "              \n",
    "\n",
    "wdg_distancing_on.observe(event_handler_distancing_on, names=\"value\")\n",
    "wdg_normal_contact_rate.observe(event_handler_normal_contact_rate, names=\"value\")\n",
    "wdg_vis_begin.observe(event_handler_vis_begin, names=\"value\")\n",
    "wdg_vis_end.observe(event_handler_vis_end, names=\"value\")\n",
    "wdg_distancing_contact_rate.observe(event_handler_distancing_contact_rate, names=\"value\")\n",
    "wdg_distancing_begin.observe(event_handler_distancing_begin, names=\"value\")\n",
    "wdg_distancing_duration.observe(event_handler_distancing_duration, names=\"value\")\n",
    "\n",
    "tabbed_graphs = widgets.Tab(children = [out_population, out_intensive_care, out_indicators, out_contact_rate, out_assumptions])\n",
    "tabbed_graphs.set_title(0, 'Population')\n",
    "tabbed_graphs.set_title(1, 'Intensive Care')\n",
    "tabbed_graphs.set_title(2, 'Indicators')\n",
    "tabbed_graphs.set_title(3, 'Contact Rate')\n",
    "tabbed_graphs.set_title(4, 'Assumptions')\n",
    "\n",
    "display(tabbed_graphs)\n",
    "\n",
    "\n",
    "control_panel = widgets.GridBox([\n",
    "    wdg_vis_begin,wdg_vis_end,\n",
    "    wdg_distancing_on,\n",
    "    wdg_normal_contact_rate,\n",
    "    wdg_distancing_contact_rate,\n",
    "    wdg_distancing_begin,\n",
    "    wdg_distancing_duration],\n",
    "    layout=widgets.Layout(grid_template_columns=\"repeat(2, 500px)\"))\n",
    "display(control_panel)\n",
    "\n",
    "\n",
    "update_graphs(wdg_distancing_on.value,\n",
    "              wdg_normal_contact_rate.value,\n",
    "              wdg_vis_begin.value,\n",
    "              wdg_vis_end.value,\n",
    "              wdg_distancing_contact_rate.value,\n",
    "              wdg_distancing_begin.value,\n",
    "              wdg_distancing_duration.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  },
  "nteract": {
   "version": "0.22.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
